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Abstract: This paper presents SPLASH, a publicly available interactive visualisation tool for 
Smoothed Particle Hydrodynamics (SPH) simulations. Visualisation of SPH data is more compli- 
cated than for grid-based codes because the data is defined on a set of irregular points and therefore 
requires a mapping procedure to a two dimensional pixel array. This means that, in practise, many 
authors simply produce particle plots which ofi^er a rather crude representation of the simulation out- 
put. Here we describe the techniques and algorithms which are utilised in SPLASH in order to provide 
the user with a fast, interactive and meaningful visualisation of one, two and three dimensional SPH 
results. 
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1 Introduction 

Sm oothed Particle Hydrodynam ics (for recent reviews 
see lMonaghan|[2005l : |Pricell2"ooi l is a Lagrangian par- 
ticle method for solving the equations of fluid dynam- 
ics. It has found widespread use in astrophysics due to 
the ability to simulate complicated three dimensional 
flow geometries and free surfaces with relative ease 
and the natural coupling with A'^— body techniques for 
self-gravitating problems. For example, SPH is used 
widely for simulations of cosmological structure forma- 
tion (e.g. Htcrik ct al. 1999; Springcl 2005), for prob- 
lems related to star (e.g . iBate et al.ll2003 :l and planet 
(e.g. iMaver et al ] |2002^ formation and in simulating 
astrophysical accreti on discs(e.g.lSmith et al.ll2007l) and 
stella r collisions fe.g. lFreitae fc Benj2005l : lDale fc Davied 
l2006l1 and pub l icly-av ailable SPH codes such as GADGET- 

2 bv ISpringell l|2005l) have found widespread applica- 
tion. 

However, visualisation of SPH data is not a straight- 
forward process, since the data is defined on a set 
of moving points which follow the fluid motion and 
derivatives are evaluated by interpolation from neigh- 
bouring points weighted by a smoothing kernel. In 
practise many authors simply present particle plots 
which are a rather crude representation of the data. 
For example the widely used and publicly available 
Tipsj{3 visualisation tool, though written for A"— body 
simulations, is often used for SPH visualisation where 
the only procedure possible is to colour the particles 
according to the value of a scalar field such as density. 

A faithful visualisation of SPH data is much more 
complicated than for grid-based fiuid codes since, for 
a smooth representation, a mapping procedure from 
the particles to a two dimensional array of pixels is re- 
quired. Using commercial visualisation packages (e.g. 
IDL) for this procedure is often inefficient because, for 
example, they require simply interpolating to a 3D grid 



first rather than mapping directly from the particles to 
the two dimensional pixel array required for a partic- 
ular visualisation. Also, given that interpolation lies 
at the heart of SPH, consistency suggests use of the 
same interpolation algorithms as part of the visuali- 
sation procedure. Because fluid particles in SPH pre- 
serve their identity, there are also certain visualisation 
procedures which are possible which cannot be used 
in an Eulerian context, such as tracing the history 
of a portion of the flow via its component particles 
and tracking of particular objects. These aspects of 
SPH visualisation give strong motivation for a dedi- 
cated software tool designed to visualise SPH data us- 
ing SPH algorithms. This paper presents the software 
design and algorithms implemented in exactly such a 
tool, which we have called "SPLASH" . 

SPLASH differs from other visualisation tools be- 
cause it is designed specifically for SPH visualisation 
and works both interactively and non-interactively (see 
the discussion relating to the software design below). 
For example IFrllO is a publicly-available tool written 
to visualise ionisation fronts in cosmological simula- 
tions (including those using particles) but allows only 
an interactive visualisation and lacks many of the fea- 
tures of SPLASH such as the ability to visualise in 
one, two and three dimensions, to select and hide par- 
ticles and to track portions of the fiow across multi- 
ple dump files. SPLASH allows plotting to both in- 
teractive and non-interactive devices allowing both a 
mouse-click driven visualisation as well as a "pipeline" 
mode for producing the same visualisation from a se- 
ries of dump files (without the need for any kind of 
scripting) . Similarly SplotclQ is a raytracing utility to 
visualise SPH simulations in a manner similar to the 
"surface rendering" technique implemented in SPLASH 
(see ij3.2p but does not allow other visualisation tech- 
niques and does not have any interactive capabilities. 



2http://home. fnal.gov/~gnedin/IFRIT/ 



-"^littp: / /www- hpcc.astro.washington.edu/tooIs/tipsy/tipsy.html ''http: / /dipastro.pd. astro. it/'^cosmo/Splotch/ 
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Other pubhcly available tools such as Supermongo 
and Gnuplot implement primitive plotting functional- 
ity at a much lower level and would require a script of 
similar length to the SPLASH source code to achieve 
similar functionality in terms of visualising SPH data 
(equivalent to SPLASH's use of the PGPLOT library 
for actually plotting the results of the rendering opera- 
tions). SPLASH can also be used to visualise remotely 
from the same location as the data is produced (e.g. 
on a remote supercomputing facility), installation on 
which is straightforward since the only requirement is 
a Fortran compiler which can also be used to compiler 
the PGPLOT libraries. Using a commercial package, 
this would not always be possible because it would 
require the remote facility to have the appropriate li- 
cense (this in particular applies to IDL). Furthermore 
many visualisation tools require some form of scripting 
to achieve the desired functionality (in IDL's case, to 
the level of an entire programming language). Since 
SPLASH is specifically tailored to visualise SPH simu- 
lations with settings changed via a series of command- 
line based menus, no scripting is required even for com- 
plicated tasks such as producing a sequence of plots 
from multiple dump files (either interactively or non- 
interactively) . 

The paper is organised as follows: In section [5] we 
discuss the basic requirements driving the software de- 
sign and present the design in detail; in ^we discuss 
the basic methods for visualising SPH data and how 
these are incorporated into SPLASH and in !j4]we dis- 
cuss the details of the interpolation algorithms imple- 
mented. Some additional features are described in ^ 
and the code's performance and memory usage are de- 
scribed in ^ A summary is given in ij?) 

2 Software design 

The basic requirements I set for an SPH visualisation 
tool (based largely on my own experience of performing 
SPH simulations) were as follows: 

1. capable of producing sufficiently annotated, ap- 
propriately labelled figures suitable for inclusion 
in research papers 

2. capable of producing a sequence of images for 
making animations 

3. capable of reading data directly from binary code 
dumps from users' SPH codes 

4. visualisation of SPH data in 1, 2 and 3 dimen- 
sions 

5. algorithms should be consistent with the basic 
SPH method 

6. should be easy to apply the same visualisation 
to different dump files (either interactively or 
non-interactively) 

7. visualisation of both scalar and vector fields de- 
fined on the particles 

8. visualisation should be interactive so the user 
can rapidly understand the data and find the 
best representation 



9. remote visualisation capability - since simula- 
tion data is often produced remotely on super- 
computing facilities from which data transfer is 
awkward and time-consuming 

10. written in a programming language familiar to 
users 

SPLASH is a program designed to meet these basic 
visualisation requirements in the most efficient manner 
possible. Each of the above requirements have strongly 
constrained the software design. For example the re- 
quirement that the visualisation be interactive means 
that simple but inefficient procedures such as inter- 
polating from the particles to a 3D grid before using 
standard grid-based visualisation techniques cannot be 
utilised. 

The basic software design which achieves all of the 
above is outlined in Figure[l] The code (written in For- 
tran 90) is built around a command-line menu struc- 
ture (designed so as to meet the requirement for remote 
visualisation) with the actual plotting performed via 
the PGPLOT graphics subroutine librarjQ (thus sat- 
isfying the requirements for production of figures for 
papers - via the postscript device drivers; for movies 
- via bitmap device drivers such as PNG and GIF; 
for interactivity - via interactive devices such as the 
X- windows driver). The use of a graphics library not 
only facilitates the easy reproduction of the same plots 
on different devices but also means that SPLASH can 
be focussed on the data input and manipulation side 
of the visualisation procedure rather than the imple- 
mentation of primitive plotting functionality. 

Plot settings are changed either non-interactively 
via a series of sub-menus accessed on the command line 
from the main menu; or interactively using the mouse 
and/or pressing particular keystrokes with the cursor 
in the plotting window (this is the "interactive mode" 
indicated in Figure [l}. 

Rather than requiring the user to convert data to 
an intermediate format (e.g. ascii files), data is read 
directly from the binary code dump files - this is a cru- 
cial requirement for rapid visualisation and makes for 
significantly reduced disk space requirements (since no 
intermediate storage is required), which can be a ma- 
jor constraint on many systems for simulations involv- 
ing > 10*' particles. The filenames are read from the 
command line, making it easy to read all files from a 
simulation by using wildcards (e.g. "splash dump*"). 
Read routines are supp lied for severa l widely used SPH 
codes (e.g. GADGET, Springej2005l : VINE.'We tzsteinI 
i2007; and Matthew Bate's SPH code, Hate 1995). Op- 
tionally, a further set of derived quantities can be cal- 
culated from the data read. For a typical SPH data 
set this would include the radius, the magnitude of all 
vector quantities and the entropy. These quantities 
appear as "extra columns" as if they had been read 
from the dump file. 

The first file listed on the command line is read on 
entry (see Figure [ij and this determines the basic pa- 
rameters used for the visualisation such as the number 
of data columns, column labels and unit settings where 
appropriate. The data is then listed by column in the 

*http: / /www. astro. caltech.edu/'^tjp/pgplot 
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Figure 1: Basic software design 



main menu, where the column number corresponds to 
the position of a variable in the data read. This means 
that any two parameters can be plotted against one 
another (for example density vs x would be plotted by 
typing 6 for the y-axis and 1 for the x-axis assuming 
column 1 contains the particle x position and column 
6 contains the density). Where two of these columns 
correspond to particle coordinates we refer to this as 
a "coordinate plot" which (provided particle masses, 
density and smoothing lengths have been read from 
the dump file) can be plotted either as particle plots 
or with a third quantity "rendered" to a pixel array. 
Thus, for example, a plot of density in a two dimen- 
sional domain is a plot of y vs x with density rendered. 
If vector quantities are present in the data (specified 
in the data read corresponding to that particular data 
format) a fourth quantity can also be plotted over the 
rendered plot in the form of an arrow plot. In 3D these 
plots can either be projection (using all particles) or 
cross section plots (using only particles contributing 
to a slice positioned in the third coordinate direction) . 
Similarly two dimensional "rendered" plots are either 
plots using all of the particles or line plots tracing an 
oblique cross section through the computational do- 
main. The interpolation procedure used to map from 
the particle data to a rendered image are described 
below and the algorithms are presented in SQl 

The plotting is directed to a particular device via 
a PGPLOT prompt. For interactive devices, the pro- 
gram then enters "interactive mode", where the user 
can manipulate the data interactively either using the 
mouse (to zoom, change colour bar limits, select and 
colour particles and move legend positions) or via keystrokes 
pressed in the plot window, giving access to a wide 
range of options such as: rotating the particles, moving 



the 3D observer, adapting plot limits, plotting smooth- 
ing circles, labelling particles, changing the colour scheme, 
adjusting the length of arrows on vector plots, set- 
ting up animation sequences, finding the gradient of 
a line and (most importantly) moving forwards and 
backwards through timesteps. For example pressing 
the space bar moves forwards to the next dump file, 
whereupon the same plot is repeated (and repeatedly 
pressing the spacebar produces a crude 'animation' de- 
pendent of course, on the speed at which data can 
be read from disk and plotted). This is indicated by 
the loop in Figure [T] which proceeds from interactive 
mode via the data read back to the "plot" step and 
finally returning to interactive mode. Where no data 
read is required the plot is simply re-plotted with the 
changed settings (perhaps recalculating the interpola- 
tion to pixels where necessary). 

A key feature facilitating the easy production of 
animations is that, when plotting is directed to a non- 
interactive device, the plotting cycles automatically 
through all of the dump files on the command line. 
This is indicated by the loop in Figure [T] proceeding 
from the "plot" step back to the data read (if the de- 
vice is non-interactive) and returning plot the same 
figure for the next dump file with settings unchanged. 

The settings for a particular plot can be saved to 
disk by pressing 's' from the main menu (see Figure [T] 
This saves a file in the current working directory con- 
taining (in Fortran 90 NAMELIST format) all of the 
current plot settings. This file is then read automati- 
cally on the next invocation of SPLASH such that plot 
settings can be restored. A "full save" (implemented 
by pressing 'S' from the main menu) saves both the 
plot settings and the current minimum and maximum 
limits set for each column (in a simple two-column ascii 
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file), so that exactly the same plot can be reproduced 
on the next invocation of SPLASH. Additional files are 
also saved where physical units have been applied to 
the data columns or animation sequences have been 
set. 

The plot settings are structured into Fortran 90 
modules which contain the parameters which may be 
changed via a particular submenu together with the 
subroutine implementing the submenu itself. Each 
settings module contains it's own namelist for those 
parameters which should be saved to disk. Thus the 
'save' operation simply saves all of the namelists in 
order into a single file. This structure means that, 
for the programmer, it is a straightforward task to 
add additional menu options affecting particular plot- 
ting functions (e.g. settings related to vector plots 
are changed in a "vector plot options" submenu and 
both the settings and the submenu are contained in the 
same Fortran 90 module. This module is then USE-d 
only in the subroutines which implement the plotting 
of vector plots, so any parameters changed via options 
in the vector submenu will be automatically available 
near where they will be used to make plotting decisions 
and automatically saved to the defaults file provided 
they have been added to the namelist). 

3 Plot types 

The "central engine" of the visualisation procedure is 
encapsulated in the "plot" step in Figure [T] An ex- 
panded outline of this step is shown in Figure[2] There 
are essentially two types of plots: particle plots or ren- 
dered plots, where a further rendered plot of vector 
arrows can be plotted on top of either of these. The 
procedure for each of these is described in turn be- 
low. Note that transformations such as log, rotation, 
3D perspective and change of co-ordinate systems are 
applied to the particle data prior to calling any inter- 
polation routines. 

3.1 Particle plots 

If rendering is not being used (ie. the plot is not a co- 
ordinate plot or no third quantity has been selected), 
the plot can simply proceed by plotting the particle 
positions directly on the plotting device (Figure [5]), 
using markers which can be chosen dependent on the 
particle types (set via the submenus accessed from the 
main menu - see Figure [T|. A simple particle plot ex- 
ample is shown in the top panel of Figure |3l Particle 
colours can be changed in a variety of ways. For ex- 
ample, selecting particles with the mouse and pressing 
keys 1-9 whilst in interactive mode colours the selected 
particles with colours corresponding to the respective 
PGPLOT colour indices. Since colours allocated to 
particles are retained in all subsequent plots, this can 
be used to select ranges of a particular parameter (e.g. 
by selecting particles on a density vs x plot) with the 
colours still appearing on a different plot (e.g. a co- 
ordinate plot of y vx x). Similarly particles can be 
coloured using data from a dump corresponding to the 
initial conditions and, provided particles retain their 




panicle data 



apply transforms 

e.g. rotate, 3D perspective 

loglO, coordinate system 
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Figure 2: Plotting pipeline 



identity between dumps, the same particles will still 
appear coloured when plotted in subsequent dumps. 

Particles can also be coloured according to the value 
of a particular quantity by setting an option which ren- 
ders via particle colours instead of interpolating to pix- 
els, although the latter method (see below) is strongly 
preferred as a method of visualisation. However there 
are circumstances where it may be desirable to see the 
actual values of a quantity on the particles themselves. 

Where the 'cross section' option has been set from 
the menu, particles are plotted in a thin slice of finite 
(although user-adjustable) thickness around the user- 
defined cross section slice position. 

3.2 Rendered plots 

"Rendered" or "pixel" plots proceed in a similar man- 
ner to particle plots but with an intermediate step 
where the particle data is interpolated to the two- 
dimensional pixel array corresponding to the viewing 
surface (Figure [2|. In 3D rendered plots are either 
projections (integration along the line of sight), cross 
section slices or surface rendered plots (see H4.3.3p . In 
2D the plot can either be a projection (a straight inter- 
polation to a 2D pixel array) or a cross section (a ID 
line plot drawn arbitrarily through the 2D domain). 
Rendered plots do not apply in the case of ID data. 
An example of a 2D rendered plot is shown in the mid- 
dle panel of Figure |3l where the same data shown in 
the particle plot (top panel) has been used. Note the 
striking difference between the visualisation using pix- 
els compared to the raw particle plot (this kind of plot 
is often used as a representation of the density field). 
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Figure 3: Top panel: A simple particle plot pro- 
duced from a 2D simulation simply by plotting all 
of the particle positions. Middle panel: The same 
data but plotted with data rendered to a pixel ar- 
ray instead of plotting particles. Bottom panel: 
with a vector plot additionally overlaid (in this 
case showing the magnetic field in the simulation) . 



One of the goals of SPLASH is to make visualisation 




Figure 4: Visualisation of a large scale star cluster 
formation calculation (Price & Bate 2007) via a 
3D rendered plot showing the density integrated 
through the z-direction. 



of SPH data in this manner a straightforward task for 
the user. 

A slight complication here is that often simulations 
contain particles of multiple types, some SPH (e.g. dif- 
ferent types of gas particle) and some non-SPH (e.g. 
sink or A'^— body particles). In this case the interpola- 
tion is performed using all of the available SPH parti- 
cles, provided plotting of that type has been turned on 
via the submenu options. Particles of non-SPH types 
can optionally be plotted on top of rendered plots (e.g. 
sink particles appear on top of a rendered plot of gas 
density) - this is indicated by the dashed pathway in 
Figure [2] An example of a three dimensional rendered 
projection (ie. showing in this case column density) 
of a large scale st ar formation calc ulation (similar to 
that described in iBate et"aD l2003h is shown in Fig- 
ure]^ where additionally a sink particle has been plot- 
ted over the rendered array. 

One further type of rendering is available in SPLASH 
for three dimensional data, which we refer to as "sur- 
face rendering" (the algorithm is described in ij4.3.3l 
below). This type of rendering provides an "optically 
thick" view of the particles (as opposed to the "op- 
tically thin" view provided by the column integrated 
rendering), showing the value of a particular parame- 
ter on the "surface of last scattering" , determined by a 
user-defined opacity which is proportional to the parti- 
cle density. Generally this type of visualisation works 
best for simulations where there is a well-defined sur- 
face and/or the range of densities in the simulation is 
not too high. An example is shown in Figure [S] show- 
ing gas temperature during the mer ger of two neutron 
stars similar to those described in iPrice fc Rosswod 
(|2006l ). The top panel shows a surface rendering near 
the start of the simulation where all of the particles 
have been used in the interpolation. The bottom panel 
shows a similar plot but where only particles below the 
midplane have been used in the calculation, giving a 
'cut-away' effect. 
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contours (instead of plotting arrows). As presently 
implemented (see ij4.2.4p this works quite well when 
the vector field is smooth but gives poor results where 
the field has strong gradients. A strongly desirable fea- 
ture for future implementation would be an algorithm 
for tracing three dimensional field lines through SPH 
particle data. 



4 Interpolation algorithms 



Figure 5: Visualisation of a neutron star merger 
calculation (jPrice fc Ros swog 2006) via a 3D sur- 
face rendered plot showing the temperature on a 
'surface of last scattering'. The top panel shows 
results near the start of the simulation using all 
of the particles, whereas in the bottom panel only 
particles below the mid-plane have been used in 
the interpolation, producing a 'cut-away' effect. 



4.1 SPH interpolation 

The heart of the SPH method (see e.g. iMonagh^nl 
Il992l : |Price| [2004: Monaghan 200^ for reviews) is the 
following identity for an arbitrary function A{r) de- 
fined on spatial coordinates r: 



Air) 



A{r')5(\r-r'\)dr', 



(1) 



where 5 is the Dirac delta function. This integral is 
approximated in SPH by replacing the delta function 
with a smooth function W with finite characteristic 
width h which reduces to a delta function in the limit 
h —> 0, giving the SPH 'integral interpolant' in the 
form 



A(r) = J A{r')W{\r-r'\,h)dr' + 0{h^), 



(2) 



where the error in the representation of A is of order 
provided the kernel function W is even and the kernel 
function is normalised such that the volume integral of 
the kernel is unity. This integral is discretised onto the 
particles by replacing the integral with a summation 
over neighbouring particles and replacing the mass el- 
ement pdr' with the neighbouring particle mass m, ie. 



3.3 Vector plots 

Vector quantities are visualised using arrow plots, al- 
though more advanced visualisations may be possi- 
ble in future. Whilst in principle an arrow could be 
plotted for each SPH particle with length proportional 
to the value of the vector on that particular particle, 
these type of plots quickly become cluttered when large 
numbers of particles are used in the simulation. Thus 
vector plots in SPLASH are implemented by first in- 
terpolating each component of the vector quantity to 
the two-dimensional pixel array corresponding to the 
viewing surface, where in 3D the plot can be an inte- 
gration of each component along the line of sight or 
where vector arrows are plotted in a cross section slice 
(depending on whether cross sections or projections 
have been selected in the menu options, also affecting 
rendered plots). 

An example of a vector plot is shown in the lower 
panel of Figure [3] where the arrows are shown overlaid 
on the rendered plot of density (otherwise identical 
to the middle panel). A preliminary feature has also 
been implemented whereby streamlines can be calcu- 
lated for the interpolated vector field and plotted as 



A(r)^^^^,W^(lr-r,|,/0. (3) 

where the subscript j refers to a quantity defined on 
particle j. The expression given above is the SPH 
'summation interpolant', forming the basis of the SPH 
approach and therefore the basis of the interpolation 
algorithms used in SPLASH for SPH visualisation. A 
normalised version of this interpolant is achieved by 
dividing the result by the interpolation of unity, given 
by 



l«^!^M^(|r-r,|,/.). 



Pi 



(4) 



Many different forms are possible for the smooth- 
ing kernel W , but the most commo nly used is the cubic 
spline kernel (see lMonaghanlll992l ): 



W{r,h) = - 



< g < 1; 

1 < g < 2; 
q>2 



(5) 



where q = \ra — Vi,\/h, v is the number of spatial di- 
mensions and the normalisation constant Oy is given 
by G\ — 2/3, CT2 ~ 10/(77r) and as — I/tt. This kernel 
satisfies the basic requirements that it is Gaussian-like 
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Figure 6: Cross section slice of density (and ve- 
locity arrows) in a neutr on star merger calculation 
( Price fc Rosswogl2006| ) showing the difference be- 
tween non-normalised (top) and normalised (bot- 
tom) interpolation. Normalised interpolation is 
turned off by default as it produces spurious ef- 
fects due to individual particles at free surfaces 
(bottom panel). 



and has smooth first derivatives which tend smoothly 
to zero as g — > 2 and is zero beyond q — 2. The 
quantity h is the smoothing length, which in most as- 
trophysical applications is a spatially variable quan- 
tity set in such a way as either to fix (either exactly 
or approximately) the number of nearest neighbours 
(|Hernauist fc Katd[l989l : iBenz et al.lll990l). or via an 



shows a cross section slice of density fro m a three di- 

mens ional neutron star merger calculation l|Price fc Rosswoe 
I2OO6I I. The top panel shows the results using a non- 
normalised interpolation whereas the bottom panel shows 
the results when the interpolated array is normalised 
(by dividing by the interpolation of unity). The nor- 
malised interpolation performs poorly at the edges, 
where the effects of individual particle smoothing spheres 
are visible. However, using a normalised interpolation 
improves the accuracy of volume rendered quantities 
on the pixels by removing effects due to the parti- 
cle distribution. Thus it is recommended that a nor- 
malised interpolation should always be used if there 
are no free surfaces. 

To avoid round-off error in interpolation calcula- 
tions (done in single precision), we write the summa- 
tion interpolant in the simpler form: 



A{v)^J2'^,AjW{r/h). 



where Wj is the dimensionless weight given by 



(6) 



(7) 



where v is the number of spatial dimensions and W 
refers to the dimensionless part of the kernel function, 
such that 



W{\r-r,\,h) = —W{r/h), 



(8) 



(ie. we have incorporated the l/ft"^ part of the usual 
kernel definition into the weight). 

With this definition a normalised interpolation is 
given by 



A{r) 



E,=i ' 



'j,AjW{r/h) 



(9) 



As an interesting aside, it is worth noting that the 
usual formula for varying the smoothing length in SPH 
codes is given by 



(10) 



where 77 is a constant and v refers to the number of 
spatial di mensions. Enforcing this relation rigourously 
(e.g. as inlSpringel fc Hernquistl2002l : |p"rice fc MonaghanI 
l2004l . [ 2007t) thus corresponds to using constant weights 
(Equation!?]) in the interpolation with the value related 
to the parameter r). Thus strictly, only knowledge of 
the (constant) weight value and the smoothing length 
is required for interpolation of any quantity in these 
codes. 



A.2 



analy tic relation to the (number) density (ISpringel fc Hernguis t 
: lMonaghanll2002l : lPrice fc Monaghanll2007i ) 



1200: 

By default the interpolations used in SPLASH are 
non-normalised. The reason for this is that, at a free 
surface the normalised interpolation (that is, using 
Eqn. [3] and dividing the result by Eqn. Q looks odd, 
whereas an interpolation using Eqn. (|3]) falls away 
smoothly. An example is shown in Figure [S] which 



Rendering of 2D data 
4.2.1 Interpolation to pixels 



Rendering of 2D data involves a straightforward ap- 
plication of Eqn. ((6)1 to the interpolation of data from 
the particles to a two dimensional grid of pixels. Thus 
we have 

A(a;,y) = ^w,^W(r/;i), (11) 
i 
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Figure 7: Interpolation of 2D data: For each par- 
ticle we perforin a loop over the pixels (in x and 
y) to which it contributes, adding the contribution 
from that particle to the pixel array. 




Figure 8: Computation of a one dimensional cross 
section through 2D data. Each particle con- 
tributes to a sequence of pixels along the section 
of the cross-section line (if any) that intersects the 
smoothing circle. 



where 



(12) 



the summation is over contributing particles and we 
take the smoothing length as 



h — max{hj, A/2), 



(13) 



that is, the maximum of the particle smoothing length 
and half of the pixel width (the latter thus being used 
generally only when few pixels are used in the interpo- 
lated plot). The interpolation is performed as a "scat- 
ter" operation from the particles, that is, for each par- 
ticle b, we find the range of pixels to which the particle 
should contribute (in both x and y) and add the contri- 
bution from particle b to all of those pixels. Note that 
this is much more efficient than attempting to perform 
the summation over particles in Equation (|lip for ev- 
ery pixel. The procedure is illustrated in Figure [7] and 
examples of 2D interpolation are shown in Figure [S] 

4.2.2 Cross sections of 2D data 

The cross-sectioning algorithm for 2D data (giving a 
ID line) is completely general and can be used for ar- 
bitrary oblique (or straight) cross sections. The cross 
section is defined by two points (xl,yl) and (x2,y2) 
through which the line should pass. These points are 
converted to give the usual equation for a line 



y = mx + c. 



(14) 



This line is then divided evenly into pixels to which the 
particles may contribute. The contributions along this 
line from the particles is computed as follows: For each 
particle, the points at which the cross section line inter- 
sects the smoothing circle are calculated (illustrated in 
Figure |8]). The smoothing circle of particle i is defined 
by the equation 



(^ 



+ (y- 



(15) 




Figure 9: Example of a one dimensional cross- 
section through 2D data, in this case showing 
the pressure distribution along a y = 0.3125 cut 
through a high resolution version of the simula- 
tion shown in Figure [31 



The x-coordinates of the points of intersection are the 
solutions to the quadratic equation 



{l+m^)x^+2{m{c—yi)—Xi)x+{xf+yf 



2cy^+c^~{2hf) 
(16) 

For particles which do not contribute to the cross sec- 
tion line, the determinant is negative. For the particles 
that do, it is then a simple matter of looping over the 
pixels which lie between the two points of intersection, 
calculating the contribution to each pixel using the ID 
SPH summation interpolant, ie. 



A(a;) = ^«;,^W(|x- 



(17) 



An example of a ID cross section through 2D data 
is shown in Figure In principle a similar method 
could be used for oblique cross sections through 3D 
data. In this case we would need to find the inter- 
section between the smoothing sphere and the cross 



= 0. 



www. publish . csiro. au /journals /pasa 



9 



section plane. However in 3D it is simpler just to ro- 
tate the particles first and then take a straight cross 
section as described above. 



4.2.3 2D vector plots 

Vector plots of 2D data are produced by interpolating 
the X— and y— components of the vector separately to 
the pixel array, which are then used to plot an array 
of arrows centred on the pixels, with length propor- 
tional to the vector magnitude. Each component is 
interpolated exactly as for scalar 2D data, ie. 

A4x,y) = ^^,A,,,W(r//i), (18) 

J 



r = Vix-x,)^ + {y-yj)\ (20) 
h = max(/ij, A/2). (21) 

The main difference between interpolation for vector 
plots and that used for rendered plots is that far fewer 
pixels are used for the arrow plots (otherwise arrows 
become indistinguishable). Thus in general the inter- 
polation for vector plots is more like a smoothing pro- 
cedure rather than an interpolation (ie. there are far 
more particles than pixels). Since we only calculate 
distances to the centres of pixel cells, this is where 
the minimum smoothing length given by l|2ip becomes 
particularly important in providing a smooth represen- 
tation of the data. An example of a 2D vector plot is 
shown in the lower panel of Figure |3l 

4.2.4 Streamlines 

For a two dimensional vector map, streamlines ( "field- 
lines") of the vector field can be plotted by integrating 
the vector field to find the stream function, contours 
of which provide the field lines. The stream function 
is given by 

^{x,y)^ J v^{x,y)dy- J Vy{x,y)dx, (22) 



such that 



9$ 

dy ' 

a* 

dx 



(23) 
(24) 



In SPLASH we compute the integral based on the in- 
terpolated velocity field on the pixel array using a sim- 
ple trapezoidal-rule integration. As presently imple- 
mented, this procedure works quite well when the vec- 
tor field is smooth but performs poorly where there 
are strong gradients present. 

4.3 Rendering of 3D data 

In three dimensions we must take either a projection 
through the whole domain or a cross section slice. 



4.3.1 Projections (line of sight integration) 

In the projection case we wish to obtain an integral 
of the rendered quantity along the line of sight. We 
begin with the 3D SPH summation interpolant in the 
form 

M^,y,z) = yZ'm.j^W{x-Xj,y-yj,z-Zj,hj) (25) 

where W is the usual (3D) cubic spline kernel ((Sjl . Tak- 
ing the integral of both sides along the line of sight 
(assumed to be along the z axis) we have 



j A(a;,j/,2)dz = ^^TTij— j W(x—Xj 



y-yj,z-zj,hj)dz. 

(26) 



This shows that the line-of-sight integration for three 
dimensions can be written as a two dimensional inter- 
polation 

A{x,y)^ [ A{x,y,z)dz = 'S2mj^Y{x-Xj,y-yj,hj). 
J j P' 

(27) 

where the 2D kernel (denoted Y) is the 3D kernel in- 
tegrated through one spatial dimension, ie. 

Y{x,y) = / W{x,y,z)dz. (28) 



For practical purposes we write Y in the form 



Y{r.y,h)^j^F{q,y) (29) 

where Qxy — r^y/h and F{qxy) is the dimensionless 2D 
kernel given by 



F{q^y) = j 



_>V(g)dq. 



(30) 



where = z/h, q^ = q^y+lz, R is the kernel radius {— 
2 for the cubic spline) and W is the usual dimensionless 
kernel function for the cubic spline, ie. 



3„2 , 3„3 



W(g) 



i(2-<?) 





+ 3^, 0<g<l; 

1 < g < 2; 
q > 2 



3 



(31) 



The integral (|30p is not (obviously) tractable analyti- 
cally (apart from at q^y = 0). However it is straight- 
forward to perform this integration numerically (for 
all qxy's from 0—^2) and store the results in a table 
for the interpolation calculation. This is the method 
adopted in SPLASH. An alternative would be to use a 
different kernel in the visualisation for which the above 
integral can be calculated analytically. 

As previously, to avoid problems with round-off er- 
ror we use the dimensionless weights defined in Equa- 
tion (O, thus writing the final interpolant (as imple- 
mented in the code) in the form 



Aix,y) 



j Ad. = Y.-,h,A,Fir.y/h), 



(32) 
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Figure 10: Computation of a two dimensional cross 
section through 3D data: Each particle contributes 
to pixels in the cross section plane that lie within 
the smoothing sphere. 



where as previously we take 

h = max{hj,A/2). (33) 

An example of a 3D column-integrated plot is shown 
in Figure showing the results of a large scale star 
cluster formation calculation (in this case showing in- 
tegrated density, ie. column density). 

In the case of vector quantities each component is 
interpolated separately in the form 

Ax{x,y) ^ J A^dz = WjhjAx,jF{rxy/hl;i4:) 

i 

Ay{x,y)^ Aydz = ^WjhjAy,jF{rxy/hl;i5) 
■' j 

where again 

ft = max(/ij , A/2). (36) 
This results in a line-of-sight integrated vector map 
which can be plotted on top of a rendered plot or as a 
standalone plot. 

4.3.2 Cross sections of 3D data 

A cross section can be taken of three dimensional data 
by summing the contributions to each pixel in the cross 
section plane from all particles within 2h of the plane 
(Figure [To]). In the implementation used in SPLASH 
the cross section is always at a fixed value of the third 
co-ordinate (ie. for xy plots the cross section is in the 
z direction). Oblique cross sections can be taken by 
rotating the particles first (the combination of settings 
can be achieved easily in SPLASH's interactive mode 
by drawing a cross section plane with the mouse, from 
which the rotation angle and slice position are auto- 
matically calculated and the cross section subsequently 
plotted). The interpolation for cross sections (e.g. in 
z) takes the form 

yl(a;,jy,2o) =^w,^>V(r//i), (37) 

j 

where 

r = ,/{x~x,Y + {y-y,Y^{zo-z,Y (38) 
h = max(/ij , A/2), (39) 



and zo refers to the position of the cross section slice. 
As previously, vector plots are achieved by interpolat- 
ing each component separately. Examples of 3D cross 
section plots are shown in Figure [6] of both scalar and 
vector fields. 

4.3.3 3D surface rendering of SPH data 

A further option for visualisation of 3D data is to use 
surface rendering (see section 13. 2p . The idea is to pro- 
duce a visualisation of the surface of a data set by per- 
forming a ray-trace through the SPH particles, with 
the density distribution giving the optical depth and 
the rendered quantity providing the colour. Thus low- 
density regions will be transparent whilst high density 
regions will be opaque. 

For a homogeneous medium the transport equation 
for a ray traced from ^ D is 

h{D) = /,(0)e-"-(^' + 5.(1 - e-^-f^'), (40) 

where is the (frequency dependent) intensity, Sv is 
the source function along the ray and r is the monochro- 
matic optical depth. The first term in (|40|l represents 
absorption (intensity decreases by e^^) whilst the sec- 
ond term represents emission. For example, at large 
optical depth (r — > oo) everything is obscured and all 
we see is the source function (ie. light emitted from D), 
whereas at low optical depth r — > the source func- 
tion contributes nothing and all we see is the previous 
intensity 1(0). 

The optical depth r is given by 




where p is the density and n is the opacity (with di- 
mensions of "cross section per unit mass"). 

For SPH visualisation the procedure is as follows. 
First of all we sort the particles in 'z' (where z rep- 
resents the distance from the observer to the particle. 
Then starting from the furthest particles, we consider 
the attenuation of a ray through each particle. Since 
what we are after is a final 2D pixel map, what we do 
in practise is take one ray for each pixel, but rather 
than taking a ray at a time (and looping over parti- 
cles), we loop over all of the particles (from back to 
front), calculating the contribution of that particle to 
all rays ('pixels') in the final pixel map. The optical 
depth through the particle is given by 

r{x,y) = j pdz, (42) 

where we have assumed that the opacity k is indepen- 
dent of z. Using the SPH summation for the density, 
we have 

r(x,j/) = K^m,Yl^(|r-r,|,ft)dz, (43) 

j 

giving just a summation involving the SPH kernel in- 
tegrated through one spatial dimension, which is the 
same as is used in the 3D projections (see 14.3.11 for 
details of how we compute this). All that remains is 
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to adjust K appropriately to give the desired surface 
position. In SPLASH an approximate value for k is 
computed according to 



o 



(my(O)d) 



(44) 



where h and rh are estimates for the average smoothing 
length and particle mass, calculated from the current 
(fixed) plot limits according toh — 0.5 * {hmin +hmax) 
(similarly for rh - the important aspect here is that 
these values do not change between dump files and 
can be restored from saved settings) and Y{0) is the 
value of the integrated kernel function ( H4.3.ip at the 
origin. The dimensionless parameter d is then a user 
defined value giving approximately the surface depth 
in terms of "number of smoothing lengths" . 

Actually, rather than computing the sum in Equa- 
tion (|43p for the whole ray, we consider just the atten- 
uation of the ray through one particle at a time, using 
the optical depth for that particle alone. Looping over 
each particle, we calculate the contribution to all rays 
(pixels) within the kernel radius 2h. That is we have, 
for each particle 

I{x,y) = /o(x,y)e-^'(^'^) + 5,(1 - e-^'(^-^)), (45) 

where 5"; is the source function (discussed below) and 
the optical depth through the particle's reach is 



n{x,y) = KmiY{x ~ x,,y - yi,h), 



(46) 



where Y is the integrated kernel function as in tj4.3.1l 
In the computation of the surface rendering, there 
are two ways of proceeding. The first option is to as- 
sign each particle an actual red, green and blue colour 
corresponding to the particle's value of the rendered 
quantity (ie. from the colour table). The source func- 
tion then consists of a red, green and blue intensity 
Si(r), Si(g), Si(^i,). Then we would add up [ie. using 
(|45|) ] the intensities in each colour (red, green and 
blue) to get final red, green and blue values at each 
pixel. The effect of this is to "blend" colours (so a red 
plus blue would make purple) , which is more like what 
happens in a real gas, but is meaningless in the sense 
that the colours produced may no longer correspond 
to those in the colour table. 

The alternative is to use a 'monochromatic' inten- 
sity - that is where the source function 5"; for each 
particle is just the value of the rendered quantity at 
the particle location. Alongside this a 'total' optical 
depth is computed along each ray. Again, we add up 
the intensities according to (|4Up , but now there is only 
a single value of I for each pixel, which corresponds to 
a final "ray-averaged" value of the rendered quantity. 
The pixel map can then be rendered in the usual man- 
ner using the ray-averaged values (which represent the 
values of the rendered quantity at the 'last scattering 
surface'). The only complication here is that we must 
make the particles optically thin to the background. 
Thus the final colours must be faded to the background 
colour (ie. black) according to the total optical depth 
computed for each pixel. The latter method is the one 




zobs - z 

Figure 11: 3D perspective: Objects at a distance d 
from the observer appear with unit magnification, 
whereas objects further away appear progressively 
smaller depending on their distance from the ob- 
server. 



used in SPLASH. However here we run into a limi- 
tation of the PGPLOT libraries, namely that the de- 
vices are limited to 256 colours, whereas we require 
256 colours also at various degrees of blackness. Thus 
at present a non-faded version is returned to the PG- 
PLOT device whilst a full (faded) version is written 
directly as a .ppm file (although without axes and an- 
notation). This is one of the limitations that would 
make it desirable to change the back-end graphics li- 
brary in future. 

An example of 3D surface rendering is shown in 
Figure [S] showing temperature in a simulation of the 
merger of two neutron stars. 



4.4 Rotation & 3D perspective 

Added perspective can be given to 3D plots by rotating 
the particles ("parallel projection") or using a depth- 
dependent 3D perspective (that is, so that objects fur- 
ther away appear smaller) . For SPH visualisation it is 
straightforward to apply these transformations to the 
particle positions prior to the interpolation procedure. 
The algorithm for 3D perspective is described below. 



4.4.1 3D perspective 

3D perspective (illustrated in Figure [TT|) is defined by 
two parameters: a distance to the observer (which we 
will call Zobs) and a distance between the observer 
and a screen placed in front of the observer (which we 
will call d). The transformation from usual x and y to 
screen x' and y' is then given by 



X = X* d/{zoBS - z) 
y = y*dl {zobs - z) 



(47) 



This means that objects at the distance d will have 
unit magnification, objects closer than the screen will 
appear larger (points diverge) and objects further away 
will appear smaller (points converge). The SPLASH 
default is a 1/10 reduction at the typical distance of 
the object (ie. observer is placed at a distance of 10 x 
object size with distance to screen of Ix object size. 
SPLASH sets this as default using the current 'z' plot 
limits as the 'object size'. 
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When using 3D perspective on interpolated plots 
the smoothing lengths of the particles are also mod- 
ified by the 3D perspective, although the smoothing 
length used to give the z length scale on integrated 
plots f Equation I32|l remains unchanged. 

5 Other useful techniques 

5.1 Fast particle plotting 

Without using hardware graphics rendering, plotting 
large numbers of particles to the screen can be quite 
slow (certainly too slow for interactive work) and pro- 
duces unnecessarily large files on vector plotting de- 
vices (e.g. postscript). Whilst one of the prime mo- 
tivations behind SPLASH is to remove the need for 
raw particle plots as a poor man's SPH visualisation, 
plots showing correlations between certain variables or 
radial profiles can still require projected plots of large 
numbers of particles. 

SPLASH uses a simple trick to speed up this kind 
of plotting by dividing the plot surface into an array 
of pixels (typically 500 x 500) and plotting up to a 
maximum of 2 particles in each cell. This results in a 
substantial speed increase with almost no loss in vis- 
ible information. Note that upon zooming the same 
criterion is applied to the zoomed-in view surface, so 
the effective resolution is increased appropriately. 

5.2 Accelerated rendering 

The slowest of the rendering techniques is the calcula- 
tion of a 3D projection through particles f ij4.3.1|l and 
the 3D surface rendering ( ^4.3.3|l since they both in- 
volve contributions from all of the particles in the sim- 
ulation, not just a subset. The former has the advan- 
tage that it can be easily parallelised (done so using 
OpenMP in SPLASH) whilst the latter is more com- 
plicated to implement in parallel (since for the surface 
rendering the contributions at each z must be added 
in order). However a simple optimisation can be ap- 
plied in both cases by taking advantage of the spherical 
symmetry of the kernel function. 

For example, considering the interpolation to the 
pixels shown in Figure [7] it is apparent that, provided 
we assume that the particle lies in the centre of the 
pixel which contains it, that the contribution to each 
quarter of the domain will be the same. Thus we can 
perform the interpolation to the top quarter of pix- 
els only and copy the result to the remaining three 
quarters, providing an in-principle speedup of 4 for 
particles contributing to large numbers of pixels. The 
caveat is the assumption that the particle lies in the 
centre of the pixel. In practise the optimisation works 
well (that is, the results are visually identical to the 
non-optimised version) except where the particles are 
regularly distributed in the domain (e.g. on a lattice 
in the initial conditions), in which case the shift in the 
particle positions can produce unwanted grid patterns 
in the interpolation. For this reason the 'accerated 
rendering' option is off by default but can be turned 
on by the user. 



6 Performance and memory us- 
age 

As discussed above, the slowest rendering techniques 
used in SPLASH are the calculation of a 3D projec- 
tion through particles and the 3D surface rendering. 
However, even these are sufficiently fast to be per- 
formed interactively. The algorithmic cost of the in- 
terpolation scales like Npart x A'^ij,, where Npix is the 
number of pixels to which each particular particle con- 
tributes. Thus larger images are more expensive. In 
SPLASH the default number of pixels is set quite low 
(ie. 200 X 200), with the idea being that a smaller 
number of pixels can be used for interactive work with 
the final step in producing the finished image to use a 
larger number of pixels. 

Whilst it is difficult to give precise timings (be- 
cause the exact time taken for the rendering depends, 
amongst other things, on how many pixels each parti- 
cle contributes to and thus how clustered the data is) , 
SPLASH is easily able to handle very large data sets 
interactively in reasonable times. For example, pro- 
ducing a rendered projection of column density from 
a three dimensional simulation containing 135 million 
particles to a 600x600 pixel image takes approximately 
55 seconds on a single processor of our local supercom- 
puter. Using the (shared-memory) parallel version on 
8 cores of the same machine takes approximately 12 
seconds. Similarly a 100 million particle simulation of 
a galactic disc takes approximately 26 seconds to ren- 
der to a 1000 x 1000 pixel image on a single processor 
and around 7 seconds on 8 cores. Using the acceler- 
ated rendering technique described above ( ^5.211 results 
in a factor of 2-3 speedup on these timings. Surface 
rendering is somewhat slower - approximately a fac- 
tor of two more expensive than a column-integrated 
projection and currently not implemented in parallel. 
However the surface rendering technique is also not as 
widely applicable to different types of simulation. 

In terms of memory use, by default SPLASH reads 
into memory an entire dump file, converted to a two- 
dimensional single precision array (where the dimen- 
sions are the number of particles x the number of 
columns). Thus for a typical "full dump file" from 
a simulation of 10'^ particles with 10 quantities (a;, y, 
z, Vx, Vy, Vz, particle mass, smoothing length, density 
and thermal energy) this would require approximately 
40Mb of storage (and hence 400Mb for 10^ particles, 
4Gb for 10* particles, etc.). Additionally a 4-byte inte- 
ger colour index is stored for each particle and tempo- 
rary memory is allocated for the two dimensional pixel 
array which is rendered to the screen, the size of which 
depends on the number of pixels chosen by the user (for 
example a 1000 x 1000 image would require a further 
4Mb). A low- memory mode for large datasets where 
memory is only allocated for those columns actually re- 
quired to make a particular plot is currently being im- 
plemented (though applicable only to binary formats 
where data columns can be read independently). In 
this mode the data is re-read from disk every time a 
different plot is made (e.g. when plotting a z—x projec- 
tion of column density instead oi an x — y projection). 
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Also plotting functionality which requires additional 
storage (such as particle colouring) will eventually be 
disabled in this mode. 

For smaller datasets, SPLASH can also be set to 
"buffer" all of the dumpfiles into memory, thus with 
memory requirements similar to the above times the 
number of dump files buffered. This provides a faster 
visualisation across multiple files for small datasets 
(since data does not have to be read from disk), pro- 
vided sufficient memory is available. 

For applications involving of the order of lO'' par- 
ticles (typical of many current SPH simulations), the 
slowest part of movie-making (ie. applying the same 
visualisation to a series of data files) is reading the data 
from disk. To speed up the visualisation in this case 
SPLASH flags whether or not each particular column 
is required for the image being produced. For data 
reads where columns can be read independently (in- 
cluding that for the GADGET code) this is then used 
to read only the required subsection of the data from 
the dump file, resulting in a much faster data read. 

7 Summary / roadmap 

In this paper we have presented SPLASH, a software 
tool for the visualisation of data from astrophysical 
Smoothed Particle Hydrodynamics simulations. The 
program is fully interactive, reads data direct from 
code dumps and can be used to visualise both scalar 
and vector SPH data in 1, 2 and 3 dimensions both to 
the screen and also to a variety of plotting devices pro- 
vided by the PGPLOT graphics library. The software 
is designed to provide the user with a rapid feel for 
the output of a simulation and a variety of efficiently- 
implemented visualisation techniques unique to SPH 
with which to represent the results. There are many 
other features of SPLASH not discussed in this paper 
which we leave the reader to discover for themselves 
(and are described in the SPLASH userguide). These 
include: 

• setting of animation sequences between frames 
in a movie 

• exact solutions to common test problems (e.g. 
hydrodynamic shock tubes, sedov blast wave) 

• transformation to different coordinate systems 
(e.g. cylindrical, spherical and toroidal coordi- 
nates) 

• particle tracking limits 

• multiple plots per page (appropriately tiled where 
so desired) 

• calculation of quantities not dumped 

Development of and improvements to the algorithms 
in SPLASH continues apace, particularly as a result of 
user feedback which has already helped to improve cer- 
tain aspects of the program substantially. In terms of 
future developments, most notable is the absence of a 
routine for plotting stream/field lines through 3D SPH 
data and it would be highly desirable to be able to do 
this efficiently directly from the particles (rather than 
highly inefficiently via interpolation to a 3D grid). 



Secondly in several places SPLASH has outgrown 
the capacities of the PGPLOT graphics library. There 
are now several other graphics libraries layered on sim- 
ilar interfaces to P GPLOT (e.g. PLPLOT and s2plot, 
iBarnes et al.|[2006l ) and migrating the back-end library 
to one of these would not represent a formidable chal- 
lenge. A more challenging alternative would be to 
move directly to OpenGL rendering, primarily for the 
speedup but also for the ease with which complicated 
3D graphics can be manipulated. However the diffi- 
culty with at least the last two of these (OpenGL and 
s2plot) at present is that inherent in the SPLASH de- 
sign is also the ability to produce, non-interactively, 
appropriately annotated graphics for use in research 
papers. Similarly the visualisation should apply as 
easily to a series of dump files (non-interactively) as 
it does to a single file (interactively or not). 

In summary, SPLASH is an efficient and capable 
software package which makes the visualisation of SPH 
data a straightforward and enjoyable task for the user. 
SPLASH is publicly available from |http : / /www . astro . ex . ac. uk/ peopl 
and is released under the terms of the Gnu General 
Public Licence. 
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